The purpose of these benchmarks is to compare several integrators for use in molecular dynamics simulation. We will use a simulation of liquid argon form the examples of NBodySimulator as test case.
using ProgressLogging using NBodySimulator, OrdinaryDiffEq, StaticArrays using Plots, DataFrames, StatsPlots function setup(t) T = 120.0 # K kb = 1.38e-23 # J/K ϵ = T * kb # J σ = 3.4e-10 # m ρ = 1374 # kg/m^3 m = 39.95 * 1.6747 * 1e-27 # kg N = 216 L = (m*N/ρ)^(1/3) R = 2.25σ v_dev = sqrt(kb * T / m) # m/s _L = L / σ _σ = 1.0 _ϵ = 1.0 _m = 1.0 _v = v_dev / sqrt(ϵ / m) _R = R / σ bodies = generate_bodies_in_cell_nodes(N, _m, _v, _L) lj_parameters = LennardJonesParameters(_ϵ, _σ, _R) pbc = CubicPeriodicBoundaryConditions(_L) lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => lj_parameters)); simulation = NBodySimulation(lj_system, (0.0, t), pbc, _ϵ/T) return simulation end
setup (generic function with 1 method)
In order to compare different integrating methods we will consider a fixed simulation time and change the timestep (or tolerances in the case of adaptive methods).
function benchmark(energyerr, rts, bytes, allocs, nt, nf, t, configs) simulation = setup(t) prob = SecondOrderODEProblem(simulation) for config in configs alg = config.alg sol, rt, b, gc, memalloc = @timed solve(prob, alg(); save_everystep=false, progress=true, progress_name="$alg", config...) result = NBodySimulator.SimulationResult(sol, simulation) ΔE = total_energy(result, t) - total_energy(result, 0) energyerr[alg] = ΔE rts[alg] = rt bytes[alg] = b allocs[alg] = memalloc nt[alg] = sol.destats.naccept nf[alg] = sol.destats.nf + sol.destats.nf2 end end function run_benchmark!(results, t, integrators, tol...; c=ones(length(integrators))) @progress "Benchmark at t=$t" for τ in zip(tol...) runtime = Dict() ΔE = Dict() nt = Dict() nf = Dict() b = Dict() allocs = Dict() cfg = config(integrators, c, τ...) GC.gc() benchmark(ΔE, runtime, b, allocs, nt, nf, t, cfg) get_tol(idx) = haskey(cfg[idx], :dt) ? cfg[idx].dt : (cfg[idx].abstol, cfg[idx].rtol) for (idx,i) in enumerate(integrators) push!(results, [string(i), runtime[i], get_tol(idx)..., abs(ΔE[i]), nt[i], nf[i], c[idx]]) end end return results end
run_benchmark! (generic function with 1 method)
We will consider symplectic integrators first
symplectic_integrators = [ VelocityVerlet, # VerletLeapfrog, PseudoVerletLeapfrog, McAte2, # CalvoSanz4, # McAte5, Yoshida6, KahanLi8, # SofSpa10 ]
5-element Array{DataType,1}:
VelocityVerlet
PseudoVerletLeapfrog
McAte2
Yoshida6
KahanLi8
Let us run a short simulation to see the cost per timestep for each method
config(integrators, c, τ) = [ (alg=a, dt=τ*cₐ) for (a,cₐ) in zip(integrators, c)] t = 35.0 τs = 1e-3 # warmup c_symplectic = ones(length(symplectic_integrators)) benchmark(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), 10., config(symplectic_integrators, c_symplectic, τs)) results = DataFrame(:integrator=>String[], :runtime=>Float64[], :τ=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]); run_benchmark!(results, t, symplectic_integrators, τs)
| integrator | runtime | τ | EnergyError | timesteps | f_evals | cost | |
|---|---|---|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | Int64 | Int64 | Float64 | |
| 1 | VelocityVerlet | 58.4098 | 0.001 | 0.00366562 | 35000 | 70002 | 1.0 |
| 2 | PseudoVerletLeapfrog | 59.2155 | 0.001 | 0.0208884 | 35000 | 105002 | 1.0 |
| 3 | McAte2 | 57.6713 | 0.001 | 0.00980374 | 35000 | 105002 | 1.0 |
| 4 | Yoshida6 | 229.896 | 0.001 | 0.0219439 | 35000 | 525002 | 1.0 |
| 5 | KahanLi8 | 530.512 | 0.001 | 0.0223923 | 35000 | 1225002 | 1.0 |
The cost of a timestep can be computed as follows
c_symplectic .= results[!, :runtime] ./ results[!, :timesteps] c_Verlet = c_symplectic[1] c_symplectic /= c_Verlet
5-element Array{Float64,1}:
1.0
1.0137927860016844
0.9873557593070432
3.935903166715239
9.082574686069826
were we have normalized the cost to the cost of a VelocityVerlet step.
Let us now benchmark the solvers for a fixed simulation time and variable timestep
t = 40.0 τs = 10 .^range(-4, -3, length=5) results = DataFrame(:integrator=>String[], :runtime=>Float64[], :τ=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]); run_benchmark!(results, t, symplectic_integrators, τs, c=c_symplectic)
| integrator | runtime | τ | EnergyError | timesteps | f_evals | cost | |
|---|---|---|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | Int64 | Int64 | Float64 | |
| 1 | VelocityVerlet | 679.942 | 0.0001 | 0.00279731 | 400000 | 800002 | 1.0 |
| 2 | PseudoVerletLeapfrog | 671.562 | 0.000101379 | 0.000301274 | 394558 | 1183676 | 1.01379 |
| 3 | McAte2 | 683.144 | 9.87356e-5 | 0.000509045 | 405123 | 1215371 | 0.987356 |
| 4 | Yoshida6 | 684.565 | 0.00039359 | 0.0128874 | 101629 | 1524437 | 3.9359 |
| 5 | KahanLi8 | 665.012 | 0.000908257 | 0.0357851 | 44041 | 1541437 | 9.08257 |
| 6 | VelocityVerlet | 378.816 | 0.000177828 | 0.00394572 | 224937 | 449876 | 1.0 |
| 7 | PseudoVerletLeapfrog | 375.186 | 0.000180281 | 0.00147033 | 221877 | 665633 | 1.01379 |
| 8 | McAte2 | 386.271 | 0.000175579 | 0.00185951 | 227818 | 683456 | 0.987356 |
| 9 | Yoshida6 | 382.206 | 0.000699914 | 0.0289537 | 57150 | 857252 | 3.9359 |
| 10 | KahanLi8 | 381.662 | 0.00161514 | 0.00651373 | 24766 | 866812 | 9.08257 |
| 11 | VelocityVerlet | 212.669 | 0.000316228 | 0.0134205 | 126492 | 252986 | 1.0 |
| 12 | PseudoVerletLeapfrog | 214.316 | 0.000320589 | 0.0113341 | 124771 | 374315 | 1.01379 |
| 13 | McAte2 | 217.365 | 0.000312229 | 0.00327997 | 128111 | 384335 | 0.987356 |
| 14 | Yoshida6 | 218.757 | 0.00124464 | 0.013221 | 32138 | 482072 | 3.9359 |
| 15 | KahanLi8 | 209.768 | 0.00287216 | 0.0680963 | 13927 | 487447 | 9.08257 |
| 16 | VelocityVerlet | 121.405 | 0.000562341 | 0.00608374 | 71132 | 142266 | 1.0 |
| 17 | PseudoVerletLeapfrog | 119.33 | 0.000570098 | 0.00520789 | 70164 | 210494 | 1.01379 |
| 18 | McAte2 | 124.652 | 0.000555231 | 0.00932556 | 72043 | 216131 | 0.987356 |
| 19 | Yoshida6 | 121.941 | 0.00221332 | 0.000567354 | 18073 | 271097 | 3.9359 |
| 20 | KahanLi8 | 119.632 | 0.00510751 | 0.181431 | 7832 | 274122 | 9.08257 |
| 21 | VelocityVerlet | 67.8505 | 0.001 | 0.0110783 | 40001 | 80004 | 1.0 |
| 22 | PseudoVerletLeapfrog | 66.9239 | 0.00101379 | 0.0224749 | 39456 | 118370 | 1.01379 |
| 23 | McAte2 | 68.7452 | 0.000987356 | 0.00814607 | 40513 | 121541 | 0.987356 |
| 24 | Yoshida6 | 69.461 | 0.0039359 | 0.0120663 | 10163 | 152447 | 3.9359 |
| 25 | KahanLi8 | 67.6691 | 0.00908257 | 0.308016 | 4405 | 154177 | 9.08257 |
The energy error as a function of runtime is given by
@df results plot(:EnergyError, :runtime, group=:integrator, xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")
Looking at the runtime as a function of timesteps, we can observe that we have a linear dependency for each method, and the slope is the previously computed cost per step.
@df results plot(:timesteps, :runtime, group=:integrator, xscale=:log10, yscale=:log10, xlabel="Number of timesteps", ylabel="Runtime (s)")
We can also consider a longer simulation time
t = 100.0 τs = 10 .^range(-4, -3, length=5) results = DataFrame(:integrator=>String[], :runtime=>Float64[], :τ=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]); #run_benchmark!(results, t, symplectic_integrators, τs, c=c_symplectic)
The energy error as a function of runtime is given by
#@df results plot(:EnergyError, :runtime, group=:integrator, # xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")
We can also look at the energy error history
function benchmark(energyerr, rts, ts, t, configs) simulation = setup(t) prob = SecondOrderODEProblem(simulation) for config in configs alg = config.alg sol, rt = @timed solve(prob, alg(); progress=true, progress_name="$alg", config...) result = NBodySimulator.SimulationResult(sol, simulation) ΔE(t) = total_energy(result, t) - total_energy(result, 0) energyerr[alg] = [ΔE(t) for t in sol.t[2:end]] rts[alg] = rt ts[alg] = sol.t[2:end] end end ΔE = Dict() rt = Dict() ts = Dict() configs = config(symplectic_integrators, c_symplectic, 2.3e-4) benchmark(ΔE, rt, ts, 45., configs) plt = plot(xlabel="Rescaled Time", ylabel="Energy error", legend=:topleft); for c in configs plot!(plt, ts[c.alg], abs.(ΔE[c.alg]), label="$(c.alg), $(rt[c.alg])s", yscale=:log10) end plt